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Abstract 

We describe our 3-dimensional numerical relativity code for the evolution of inho- 
mogeneous cosmologies. During the evolution the constraint equations are monitored 
but not solved. The code has been tested against perturbation theory with good re- 
sults. We present some runs of inhomogeneous inflation with strongly inhomogeneous 
initial data. 
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I. INTRODUCTION 

Inflation[l,2] has become the favorite explanation for both the large-scale homogeneity 
and the small-scale inhomogeneity of the universe. The recent observation of the cosmic 
microwave anisotropy by the COBE satellite[3] is consistent with the inflationary universe 
scenario. Most work on inflation has been in the context of either a flat background 
spacetime or perturbation theory. A consistent general-relativistic treatment of strongly 
inhomogeneous inflation appears to require numerical relativity. 

Previous numerical relativity calculations of inhomogeneous inflation have been 1- 
dimensional[4-7]. We have written a 3-dimensional numerical relativity code for this pur- 
pose. The code consists of two parts. The first part solves the initial value problem. 
We described it in Ref. [8], where the motivation and background of this project are also 
discussed. 

Here we present the second part, the evolution code. The main focus of this paper 
is on the code rather than physical conclusions about inhomogeneous inflation. Thus we 
discuss the evolution equations (Sec. II), numerical techniques (Sec. Ill), and code tests 
(Sec. IV). In Sec. V we present a couple of strongly inhomogeneous inflation runs. 

II. EVOLUTION EQUATIONS 

We derive evolution equations for general relativity with scalar field dynamics in a 
form suitable for numerical solution. They will be in the 3+1 formalism [9] and for a 
general gauge. 

The matter source is a scalar field. The covariant form of the field equation is 

<7^ ;M „ = V'(0), (2.1) 

where V(<p) is a potential function, the form of which we are not restricting. The prime 
denotes derivative with respect to the argument. Also, the semicolon denotes covariant 
derivative with respect to the 4-metric = h^ — n^n^, and _D M is the covariant derivative 
with respect to the 3-metric h^; n M is the unit normal vector to the 3-space slice, n^n^ = 
— 1, and n M /i Mi , = (Greek letter indices range and sum over 0-3; latin letter indices over 
1-3). Since 

h^4>,^ = D^D^ + Kn^(f> (2.2) 

and 

= n v d v r) + n v d^{K» + a»n u ) (23) 

= n v d v r) - -Wductdvfa 
a 

Eq. (2.1) becomes 

D^(j) + K V - n v d u r] + hi^d„ad v (j) = V' (</>). (2.4) 
Here we have defined a new variable 

7/ = (2.5) 
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to replace a second-order equation (2.1) with two equations (2.4,5) first order in time. 
Hereafter we use 3+1 coordinates (i, x 1 , x 2 , x s ) = (i, x, y, z), with the metric[9] 

ds 2 = -a 2 dt 2 + h ij (dx i + p i dt)(dx j + P J dt). (2.6) 

The scalar field equations just derived are now 

-(d t - p k d k )rj = -=di{Vhh i3 ' drf) + K v + -h^d.ad^ - V'(<f>) (2.7) 
a ,JJi a 

and 

-(d t - p k d k )<t> = V - (2.8) 
a 

The 3+1 form of the Einstein equations is 

R + K 2 - K tJ K lJ = 16nGp H (2.9) 
DiK) - DjK = (2.10) 

d t K) = -D l D 3 a + a[R) + KK) - ZttGS) + AnG5){S - p H )] 

+ p k d k K) - K)d k p + dj[3 k Kl (2.11) 
dthij = -2aK i:i + Di(3j + Djfy. (2.12) 

Here R and Rj are the Ricci scalar and tensor of the 3-metric, and G is Newton's constant. 
Eqs. (2.9) and (2.10) are the Hamiltonian and momentum constraints. We do not concern 
ourselves with the solution of the constraint equations in this paper, as we have discussed 
it elsewhere [8]. 

For a scalar field, the source terms in Eqs. (2.9-12) will be 

p H = IWdiWrf + \i] 2 + V((f>) (2.13) 

= -rjh ij d j( j) (2.14) 

5« = h ij [-\h kl d k <j>di<l> + \r} 2 - V ((/>)] + h ik h3 l d k (f>di(f> (2.15) 

S = SI = -\h^ d l( j)d 3 (j) + \r} 2 - 3V(<f>). (2.16) 

The evolution equation for the extrinsic curvature, Eq. (2.11), becomes 



-(dt - l3 k d k )K) = KK) - -h ik (djd k a - F l jk dia) + R) 
a a j j 



8ttG \h ik d k (j)dj(j) + S)V{4>)] + ^(^/? fe ~ K k d k (3 l ) 



(2.17) 



We shall evolve the trace K = K l - and the traceless part A 1 - = K) — ^5jK separately. The 
evolution equation for the trace is 

-(d t -(3 k d k )K = K 2 -^ 7 =d l (Vhh^d J a)+R 

a aVh (2.18) 

- ZirGWd^djcj) - 24irGV{c/)). 



Using 



Difij = dif3j - \f3\djhn + dihij - tyk^), 



the 3-metric evolution equation, Eq. (2.12), becomes 

hfr - P k d k )h tJ = -2Kij + ^(hkidjP* + h kj di(3 k ). 



(2.19) 



(2.20) 



It turns out to be useful to have a separate evolution equation for the square-root of the 
determinant of the 3-metric, y/h. Using 



with Eq. (2.12) we get 



yh 



-(d t - /3 k d k )Vh = -VhK + -Vhd k (3 k . 

a a 



(2.21) 



(2.22) 



To have a complete set of evolution equations we need to specify equations for the 
lapse a and the shift i.e., to choose a gauge. The simplest choice is the synchronous 
gauge, 

a = 1 

(2.23) 

/3 l = 0. 

Another choice particularly straightforward to implement in this context is the harmonic 

gauge, 



a 



id t -(3 k d k )a = -aK 



-(d t -p k d k )F = -a 
a 



h ij ^ + -=di(Vhh* j ) 
« Vh 



(2.24) 
(2.25) 



The evolution equations, Eqs. (2.7,8,17,18,20,22), are valid for any gauge and we will 
retain this property as we code them. We do not attempt to simplify them by choosing a 
particular gauge. 

We code these evolution equations in conservative form. The conservative form of 
Eq. (2.22) is 

d t Vh - d k ((3 k Vh) = -aVhK. (2.26) 



From this follows a general conversion rule: If the original form is 

(d t -p k d k )f = g, 



(2.27) 



then the conservative form is 



d t f-d k (p k f) = g-aKf, 
4 



(2.28) 



where 

f = y/hf, (2.29) 



etc. 



Our final set of equations is 

d t K - d k (f3 k K) = -di(h ij dja) + aVhR - 8TrGah ij drfdjcj) - 24nGaVhV((PX2.30) 
d t A) - d k ([3 k A)) = -~h ik (djd k a - T l jk dia) + ah ik R kj - %irGah ik d k <t>d j <l> 

~ I 6 j [-d k {h kl dia) + aVhR - 87rGa~h ki 'd^drf] 

+ d j p k A{ - dkFA) (2.31) 

d t fj - d k {[3 k fj) = di(ah ij d j( /)) - aVhV'(cf)) (2.32) 

d t Vh - d k {[5 k Vh) = -aVhK (2.33) 

dthij - d^Kj) = -%aK~h l3 - 2aA l J h ll + dj0 l hu + d^kj (2.34) 

d t 4> - d k (f3 k 4>) = -aK4> + afj. (2.35) 

The harmonic gauge equations in conservative form are 

d t u! - d k ((3 k u) = (2.36) 
dti? - d k ((3 k (3 l ) = -aK(3 l - a 2 



a 



(2.37) 



where u = 1/a, or ui = y/h/a. 

We need initial data that satisfies the constraint equations (2.9) and (2.10), which we 
can write as 

R + \K 2 - A)A\ = STiG[h 13 d l (f)d j (t) + r] 2 + 2V ((/))] (2.38) 
diA) + Y l u A l 3 - F^A} - \djK = -8nGf]d j( p. (2.39) 

We discussed the procedure for solving the initial value problem in Ref. [8] . In addition to 
these numerically solved initial data, we have also used data that satisfies the constraints 
trivially. 

This latter method, previously used by Goldwirth and Piran[6], involves two scalar 
fields to achieve a momentarily flat space for the initial slice. One of the scalar fields is the 
inflaton field, 0, which can be set arbitrarily according to the problem we wish to solve. 
The other field, 0^, is an auxiliary field, with V(4>r) = 0. We want the right hand side 
of Eq. (2.38) to be constant, and that of Eq. (2.39) to be zero on the initial slice. This is 
achieved by setting 

rj = 

0i? = (2.40) 
m = [2 PH - WdiWrf - 2V (</>)] 1/2 , 
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where pu = const. Equations (2.38) and (2.39) can now be satisfied with 



hij = 8). 



(2.41) 



The initial slice will thus have a homogeneous total energy density, but with an inhomo- 
geneous composition. The energy density becomes inhomogeneous immediately when it is 
evolved. In a typical run, 0# soon redshifts away whereas the vacuum energy V(4>) keeps 
the inflaton field dynamically important. 



III. NUMERICAL TECHNIQUES 

Replacing the continuous spacetime with a discrete grid requires a decision about the 
relative positioning of variables on the grid. As discussed below, we have divided the 
variables into two sets staggered in time. This division arises naturally from the form 
of the equations, although some terms in the equations violate this structure and require 
interpolation. 

In the previous section, Eqs. (2.30-35) are written in the order they are actually 
evolved in our code. The lapse is calculated before Vh, and the shift before K. This 
ordering is determined by what quantities are needed to evolve which quantity. The quan- 
tities a, y/h, hij, and 4> are evaluated after the full time-step, i.e., at times 

t n - 1 ,t n = t n - 1 + At n ,..., 

where At n is the n th time-step. The quantities (3 % , K, A* , and fj are evaluated at half-time, 
i.e., at times 

t n ~ l + \At n ,t n + \At n+1 , . . . . 

With the possible exception of lapse and shift, this scheme does not require any extrapo- 
lation in time. All quantities that are needed to evolve a particular quantity have already 
been evolved to the time level at which they are needed. The quantities a, f3 l , \fh, and 
K are needed at both full- and half-time, so they are interpolated when needed. Depend- 
ing on the choice of gauge, the calculation of a and f3 l at a new time level may require 
extrapolation of some quantities in time, and should then be iterated. 

It is a common practice to have a similar staggering in space. A gradient of a variable 
would naturally lie between the grid points where the variable is located. In a 3-d code 
this becomes complicated. Picture the grid dividing the space into small cubes, or zones. 
Some variables might naturally lie at zone centers, others at zone faces, and yet others at 
zone edges. The grid will have 3 times as many faces and edges as zones (in 3 different 
orientations). Typically the different components of a 3 vector would lie on differently 
oriented faces. If the equations respect such a structure, the extra coding effort will pay 
off by giving a good accuracy with a small number of zones. We, however, did not find 
our equations to naturally fit into such a scheme. Because so many terms would have to 
be interpolated, it appeared unlikely that this would lead to much improved performance. 



6 



Therefore we decided to position all the variables the same way in space. We refer to these 
locations simply as grid points. To obtain a value of a spatial derivative at a grid point we 
thus need to difference over two grid spacings. 

For simplicity and generality, we are using Cartesian coordinates. We have coded for 
nonuniform grid spacings in x, y, and z. However, to maintain second-order accuracy in 
space with a nonuniform grid the spacing should not vary too steeply. A variable will lie 
at a gridpoint (x, y, z), where, e.g., x takes on any of the values x\, . . . ,xn x , and 

x n = x n - 1 + Ax n . (3.1) 

The first and last grid points in each direction, e.g., at x\ and xn x , are dummies. Vari- 
ables are not evolved at these grid points but are set after the evolution according to the 
boundary condition. Thus the 'physical' grid has only (N x — 2)(N y — 2)(N Z — 2) points. 
Periodic boundary conditions are 

f(x Nx ,y,z) = f(x 2 ,y,z), 
f{x x ,y,z) = f(x Nx - U y,z), 

etc. Reflective boundary conditions for tensors are different than for scalars, because some 
components change signs, e.g., 

P 1 (x 1 ,y,z) = -(3 l {x 2 ,y,z), 
h 12 (x 1 , y, z) = -h 12 (x 2 , y, z) 

The largest piece of the evolution code is the calculation of the Ricci tensor R^. It 
is calculated from the 3-metric h^ and is needed for evolving the extrinsic curvature. In 
the process we obtain the connection coefficients T l jk , also needed for evolving A 1 -. (r* fc is 
also needed for checking the momentum constraint, and R = R\ is needed for checking the 
Hamiltonian constraint). These are computed directly from their definitions 

Y) k = \U\d k h l3 + d 3 h lk - dih jk ) (3.4) 

and 

r> — Q -pk o -pfc i -pfc -pZ -pfc -pZ 
-fly = Cfcl — OjL ik -+- 1 lk L ij — 1 ijL ik 

= \d k h kl (djhH + dihij - dihij) 
+ ±h kl (d k djhii + d k dihij - d k dih i3 ) 

- \d 3 h kl d % h lk - ^djdthtk + rt k i% - v%v\ k 

(We first invert hij to obtain /V- 7 ) . We use the latter form for Rij to avoid taking derivatives 
of r* fc , which already involves derivatives. With our grid structure we obtain second 
derivatives more compactly by calculating them directly. 

There are 6 components of Rij and 18 components of T l jk to be calculated, with double 
sums. This could have been coded with loops over the indices, relying on the compiler to 
unroll the summation loops to permit vectorizing. However, there are complications due to 
the symmetries in indices. Also, many of the terms cancel out with certain combinations 
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(3.2) 



(3.3) 



(3.5) 



of indices. Therefore we have coded in all components separately, with the sums written 
out. We actually wrote a small program to produce this long (almost 600 lines of Fortran) 
but repetitive piece of code. 

The evolution equations are mostly rather straightforward to code. Since we are 
evolving the 'twiddled' quantities — see Eq. (2.29) — which are often needed 'untwiddled', 
we have to be careful to multiply and divide by \/h when needed. To make this possible, 
\fh has to be evolved to a new full-time level first. 

The present version of the code does not yet include the advection terms dk((3 k /), and 
is therefore restricted to j3 l = 0. Also, the right-hand-side terms involving j3 l have been 
left out in the A 1 - and hij equations. This allows us to evolve A 1 - without matrix inversion. 

The evolution of involves matrix inversion since the right-hand side of Eq. (2.34) 
contains other components of this tensor in the — 2a A l -hn term. With (3 % = 0, the dis- 
cretized form of this equation (taking us from t n to t n+1 ) can be written as 

H n + 1 = H n -±aAtM(H n + H n+1 ), (3.6) 

or 

H n+1 = (1 + ±aAtM)-\l - \aAtM)H n , (3.7) 



where H = (/in, h\2, ^13, ^22, ^23, ^33) is a 6-component vector, and M is a 6 x 6 matrix, 
whose components are linear combinations of the extrinsic curvature components. A matrix 
inversion at each grid point would be rather time consuming. As we only need second-order 
accuracy in At, we replace Eq. (3.7) by 

H n+i = (1 _ a AtM + \a 2 At 2 M 2 )H n , (3.8) 

which requires mere matrix multiplication, and does not prevent vectorization over the 
grid. 

A similar procedure will have to be utilized in the evolution of A % - once nonzero shift 
is included. Presently we evolve each component separately. Although only 5 of the 
components are independent, we evolve all 9 of them since reconstructing the remaining 
components from the minimum set of 5 is nontrivial. (Since A\ = 0, we could easily skip 
evolving one of the diagonal components) . The covariant components are symmetric, so if 
we rewrote Eq. (2.31) to evolve A^j instead of A 1 -, there would be only 6 components to 

evolve. But this gives a nonlinear equation, whereas Eq. (2.31) is linear in A 1 -. 

We are presently running free evolution, i.e., the constraint equations are not solved 
after the initial slice. To see how accurately the constraints are maintained, we compare 
the left- and right-hand sides of Eqs. (2.38) and (2.39) during the evolution. 

We use three conditions to determine the time step. The first is the Courant condition. 
The code requires three time steps to carry information from the grid point (xi,y m , z n ) 
diagonally to z n +\). The time step must be short enough that this grid velocity 

is not less than the speed of light. To satisfy this condition everywhere, we need to find 
the shortest physical distance As, where 

As 2 = Ax* Ax*, (3.9) 
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between neighboring points. To be more precise, we set 

At c = Cmin|^), (3.10) 
\aa J 

where d = 1,2, or 3 depending on whether the neighboring point differs in one, two, or 
three coordinates, and C < 1. In practice we have used C = 0.7. 

The second condition is an expansion condition. According to Eq. (2.33), the local 
volume expansion rate is given by —aK. To keep expansion per time step below a certain 
amount E, we set 

At E = E/ max(a\K\). (3.11) 

The third condition is a smoothness condition. To prevent the time step from increasing 
too rapidly, we set 

At S = (1 + S)Atprevious. (3.12) 

The time step is then the smallest of these three, 

At = min(At c , At E , At s ). (3.13) 

The Courant condition gives a shorter time step when the number of grid points is 
increased. In the runs presented here, we adjusted the other two conditions with grid size 
so that halving the grid spacing also halved the time step consistently over the whole run. 
Thus we have used E = 0.1(32/iV) and S = 0.01(32/JV), where N is the number of grid 
points in one direction. 

In practice, the controlling condition was the Courant condition in the early part of the 
run, the expansion condition during most of the inflation, and the smoothness condition 
after inflation. The time-scale of the post-inflation oscillations of (p is not related to any of 
these three conditions. Indeed, when we ran long enough, the time step became too long 
to resolve these oscillations. We experimented with a fourth condition to control this, but 
could not easily find one that is both general and simple. Instead, we just chose a small 
enough value for S to resolve the first oscillations, as we were not interested in the later 
ones. In running different problems it may be useful to experiment with the time-step 
control, to optimize it for each case. 

IV. CODE TESTS 

A. Homogeneous tests 

The simplest test is a homogeneous one. We did a homogeneous chaotic inflation[10] 
run with the potential V((f>) = |A<^ 4 , A = 10~ 3 , and the inflaton field initially at = 5mp\. 
The perturbation runs below are perturbations around this homogeneous run. Figures 1 
and 2 show the time evolution of the inflaton <fi and the scale factor a, respectively. The 
code deviates from the exact results by less than 10~ 3 . Another run with a different time 
step showed the error to be quadratic in time step. We also ran a Kasner model, where 
the expansion rates along the three axes were all different, and obtained similar accuracy. 
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B. Perturbation tests 

The perturbation theory for inflation has been discussed by Brandenberger et al. 
[11]. The unperturbed spacetime is assumed to be spatially flat. We shall work in the 
synchronous gauge. The metric and the inflaton field can be written as 

ds 2 = -dt 2 + a(t) 2 |[l + A(t,x)]6ij + B(t,x) jij ^dx i dx j (4.1) 

and 

0(£,x) = o (t) + <ty(*,x), (4.2) 

where 5(f), A, and B are small perturbations. The th order, or homogeneous, equations 
are 

# 2 -(£) 2 = ^[^ + n<M], (4.3) 

2^+(^) 2 = -8nG[^ 2 -V(Ml (4.4) 

and 

00 + 3^00 + ^(00) =0, (4.5) 

where Eq. (4.3) is the Hamiltonian constraint. The momentum constraint vanishes identi- 
cally to th order. 

We expand the perturbations in plane waves 

50(t,x) = ^50 fc (t)e lk - x , (4.6) 

k 

etc., where k is the comoving wave vector. The first order perturbation equations for the 
modes 8<p k {t), A k (t), and B k (t) are then 

3A k - k 2 B k + 3H(3A k - k 2 B k ) + k 2 a~ 2 A k = -24nG[<j> 5<p k - V'Mdfa] (4.7) 

A k + 3HA k = -8nG[<p 5<p k - V'i^Hk] (4.8) 

-k 2 a~ 2 A k - H(3A k - k 2 B k ) = -8nG [<j> 5<j> k + ^'(0o)^fe] (4.9) 

A k = SirGfoWk (4.10) 
5'<p k + 3H5<p k + V"((t)o)6(t) k = -k 2 a~ 2 (j) k - ±M3A k - k 2 B k ). (4.11) 

To obtain 'exact' perturbation theory results we solved these ordinary differential 
equations with the Runge-Kutta method. The metric evolution equations (4.4), (4.7), and 
(4.8) can here be ignored, and the metric (a, A k , B k ) can be evolved using the constraint 
equations (4.3), (4.9), and (4.10). 

For the perturbation runs with the code we need the relation between the code vari- 
ables K and Aj, and the perturbation theory variables A and B. From the metric evolution 
equation (2.12) we obtain 

K = -3H, (4.12) 
5K = |V 2 S, (4.13) 

Aj = -iS,y + i<SyV 2 S, (4.14) 
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with K = K + 5K. The scalar 3-curvature is 

R = -2a" 2 V 2 A 



(4.15) 



The initial data has to satisfy the perturbation theory constraint equations (4.9) and 
(4.10). This is accomplished by choosing arbitrary values for 0o, <fio, S(/>, S(/>, A, and B, and 
solving for A and B from Eqs. (4.9) and (4.10). For our test runs we have chosen initial 
data of the form 

a = 1, 
4>q = 5mpi, 
00 = 0, 



ikx 



(4.16) 



5cj) = 0, 
A = 0, 
B = 0. 

For this case, the constraint equations (4.9) and (4.10) give 

A = o, 

8ttG t// . 



Bh = 



Hk 2 



(4.17) 



The initial values for the extrinsic curvature variables are thus 

K = -3tf , 



SK 

A\ 



(4.18) 



where 



if 



8tvG 



-,1/2 



5(p k e 



ikx 



(4.19) 



The particular initial data we have used is a superposition of three perpendicular 
plane waves, 

<ty(*o, x) = <tyi(t )e ifclK + <ty 2 (*o)e i * 31 ' + 5<p 3 (to)e lk3Z , 



with small and equal amplitudes 

5Mto) = 6<h(t ) = 5Mto) = 10" 6 , 
but different wavelengths Li = 2n /ki, 

L x = 0.2H~\L 2 = 0.3H- 1 , L 3 = OAH' 1 . 
11 



(4.20) 



(4.21) 



(4.22) 



Here the subscripts i = 1,2,3 denote the three different modes, rather than components 
of a vector. Our grid had one wavelength in each direction. To get the same number 
N x = N y = N z of grid points per each wavelength, we have used different grid spacings 
Ax 1 in each direction. 

We studied the convergence of the code results towards the perturbation theory results 
as the number of grid points is increased (the time step and grid spacing is decreased). 
The results are shown in Figs. 3, 4, and 5. These show the evolution of the amplitude of 
the three perpendicular perturbation waves in <f>, A, and B. 

Initially, the waves are well inside the horizon. (By 'horizon' we refer to the Hubble 
length, as is common among astroparticle physicists). The Hubble constant at the initial 
time is 

*-(£r(£) , -" ta * (423 » 

The initial values of ki/H were 10-7T, 207r/3, and 5n. We have arbitrarily set the initial 
time as to = £pi- As the universe inflates, the waves exit the horizon one by one. Defining 
the time of exit by ki = H, this happens at t = 2.20£pi, 2.33tpi, and 2.51tpi. 

The behavior of the perturbations change markedly as they exit the horizon. Inside 
the horizon, the initial <j> perturbation oscillates and these oscillations are damped by the 
expansion, all three wavelengths at the same rate. There was no initial perturbation A, 
B, in the metric, but these are now induced by the <fi perturbation. The A perturbation 
oscillates with a growing amplitude, a quarter phase behind 5<p. The B perturbation ex- 
hibits monotonic growth with a small oscillation superimposed. The metric perturbations 
don't have time to reach the cf> perturbation amplitude, before the perturbations exit the 
horizon. 

Outside the horizon the oscillations cease. The (f> perturbation decreases while A grows 
and B stays constant. At the end of inflation the growth in A ceases, whereas 8(f) shows 
anharmonic oscillation with </> at the bottom of V ((/)). 

The discretization error exhibits a quadratic behaviour when we increase the number 
of grid points from 16 3 to 48 3 . The errors in B are small, as are the errors in A and 
5(f) while inside horizon. The larger errors in 8(f) and A outside the horizon are due to 
the switch from oscillating to monotonic behaviour, when a small phase error leads to a 
relatively large amplitude error. Even then, the errors in the final values in the 48 3 runs 
are less than 3%, after a volume expansion by more than 10 100 . 



C. Memory and time requirements 

The code has 56 variables per grid point. These are \/h, hij (6), /i lJ (6), T l - k (18), 
Rij (6), .R, a, f3 l (3), K, A 1 - (9), 0, rj, and t)r. In addition, 30 memory locations per 
grid point are used for work space or auxiliary variables. The code is not optimized for 
memory usage but rather for simplicity, so memory savings could be achieved at the cost 
of making the code more complicated, less readable and more time-consuming. Thus the 
code requires 86 words of memory per grid point, or about 24 Mwords for a 64 3 run. In 
practice, the largest grid that fit in the 32M queue on the Cray-2 at NCSA was 70 3 . 

The code runs at 140 Mflops on a single processor, and takes about 0.02 ms of com- 
putation time per grid point and cycle. Thus a 64 3 grid runs at about 700 cycles per hour. 
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To run the 70 e-foldings of inflation required to obtain the observable universe from one 
pre-inflation Hubble volume takes almost 10 hours of Cray-time with this resolution. The 
most time is spent on calculating the connection coefficients T l - k and the Ricci tensor R^j, 
about 43% of total time. The second most time is used in the evolution of the extrinsic 
curvature tensor A*, 11%. 

V. INHOMOGENEOUS INFLATION RUNS 

We now present some runs of inhomogeneous inflation. These runs use the flat initial 
data described in section II. A run using initial data obtained with the initial data solver 
was presented in Ref. [12]. 

The requirement that the initial time slice has homogeneous total energy density (as 
implied by flat initial data) gives an upper limit for the variation of the inflaton field (f> 
within a Hubble length. From 

i(V</>) 2 <p=^tf 2 (5.1) 

we get 

Thus inhomogeneities that begin inside the horizon have small amplitudes. Runs where 
we had large initial variation in the inflaton field have initial grid lengths L equal to many 
Hubble lengths. 

These runs are of chaotic inflation. The potential is V(4>) = |A0 4 , with A = 10~ 6 . 
We chose initial data of the form 

2 x 

4>(t ,x) = 4> + 8(j) ^2 j^smxismy m smz n , (5.3) 

l,m,n=l 

where x\ = 2nlx/L + 9 x i, etc., with Q x \ random phases. The initial Hubble constant was 
chosen to be H = O.lmpi, and </>o = 5mp\. 

We present a small-scale run with L = H -1 , 5<p = 0.0125mpi, and a large-scale run 
with L = 32H~ 1 , dtp = 0.4mpi. The phases for these particular runs were Q x \ = 3.83, 
9 x2 = 2.59, yl = 3.02, 6 V 2 = 1-83, zl = 3.19, and 6 z2 = 1.33. The same phases were used 
in both runs. In the initial data, the minimum <p value turned out to be <po — 3.350, and 
the maximum (po + 2.15(p. The maximum value of V</> was within a factor of two from the 
upper limit. In Fig. 6 we show the regions of small and large <p. 

The results are shown as 3-dimensional contour plots of the scalar variables 0, Vh, 
K, and R. These are supplemented by 1-dimensional plots showing the quantities along 
a reference line at different time slices. This reference line runs parallel to the y-axis and 
goes through the initial (p minimum. Figure 6 applies to the initial data of both the small- 
and the large-scale runs, although the amplitude of these 4> variations were of different 
magnitude. 

We discuss the small-scale run first. The grid length was initially set equal to one 
Hubble length. Thus the inhomogeneities are at first well inside the horizon. The dynami- 
cal time scale of the inflaton field is therefore shorter than the expansion time scale, given 
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by the Hubble time. At the initial slice the energy density is dominated by the kinetic 
energy of the <j) an d 4>R fields, accounting for 84% of total energy. 13 % of the energy is 
in V(4>) and 3% in field gradients (although at the point of maximum V0 they account 
for more than half of the energy density). The initial Hubble time is tn = lOtpi. We 
arbitrarily set the initial time as to = t#(to)/2 = 5£pi. The time step is at first controlled 
by the Courant condition, and thus we get several hundred time steps per Hubble time. 
We show a number of contour plots of the slice at cycle 256, or t = 9.2tpi. Of the total 
energy in this slice, 52% is kinetic, 46% potential, and 2% gradient. The Hubble time is 
now roughly 18tpi, so this slice is about a quarter of a Hubble time from the beginning. 

As the inhomogeneities evolve, they oscillate and their amplitude is damped, see Fig. 7. 
After a quarter of a Hubble time of evolution, the regions of small and large <fi are quite 
different from what they were initially, see Fig. 8. By the time the inhomogeneities have 
exited the horizon, the inflaton field has become much more homogeneous than it was 
initially. 

Since we started with flat initial data, y/h, K, and R were all homogeneous at the 
initial time (and R was zero). The inhomogeneities in induce inhomogeneities in these 
curvature quantities. 

We show regions of small and large expansion \fh after a quarter of a Hubble time 
in Fig. 9. Regions of both large and small initial <p end up as regions of small expansion. 
The region of largest expansion forms a shell surrounding the initial <p minimum. This is 
partly due to the large gradient energies in this region at early times, and partly due to the 
inflaton field happening to have large potential values in this region during the time when 
most of the expansion inhomogeneity is generated. After a volume expansion of about 100, 
further expansion is rather homogeneous, and the y/h profile stays the same, see Fig. 10. 

The extrinsic curvature K is related to the rate of change in \fh. In Fig. 11 large abso- 
lute values of K in the region of large \fh are seen. Later K becomes rather homogeneous 
(see Fig. 12). Regions of largest intrinsic curvature R (Fig. 13) seem to be at or around 
regions of small expansion. This is partly because in this run those are more pronounced 
than regions of large expansion, partly because expansion decreases R by increasing the 
curvature radius. Figure 14 shows R along our reference line at different times. 

We show only the early part of this run, since the later evolution is less interesting. 
The inflaton rolls down its potential staying rather homogeneous. Inflation ends as 4> 
reaches the bottom and begins to oscillate. 

In the large-scale run the inhomogeneities are outside the horizon. They do not 
oscillate. Rather, they maintain their shape without damping. After a few Hubble times 
the regions of low (Fig. 15) and high (p look very much the same as initially (Fig. 6). 
Figure 16 shows 4> along the reference line at various times through the entire run past 
the end of inflation. We see that aside from the flattening of the sharp initial minimum, 
the inhomogeneity maintains its shape. Because of the larger scale, the Courant condition 
allows a larger initial time step for this run. 

The spatial pattern of the other scalar quantities follows that of <p. Larger values of <fi 
lead to faster expansion. At first the effect of the other scalar field 4>R (see end of section 
II) shows, but later the regions of the most expansion match those of large <fi rather closely 
(Fig. 17). Since the regions of large <p and high inflation stay the same, the inhomogeneity 
in the expansion becomes steeper and steeper (Fig. 18). When inflation ends, some regions 
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have expanded by a volume factor of more than 10 130 , while others have expanded by less 
than 10 75 . 

For the extrinsic curvature we show just the 1-d plot (Fig. 19). The 3-d plots would 
look intermediate between the <f> and y/h plots. Note the correspondence between Figs. 16, 
18, and 19. Larger 0, and thus V((f>), means faster expansion, i.e., larger (more negative) 
K. This then builds up as a larger cumulative expansion \fh. 

In Fig. 20 we show regions of large positive and negative 3-curvature R at cycle 64, 
after half a Hubble time of evolution. The largest curvature is at the region where the 
expansion has been the smallest. This region of large positive curvature is surrounded by 
a shell of negative curvature. After more inflation, the curvature decreases but the relative 
contrast between large- and small-curvature regions increases, so that later contour plots 
would pick up only the largest-curvature region near the back corner. 

The Hamiltonian and momentum constraints were monitored during the evolution. 
During the small-scale run the maximum errors (root-mean-square relative difference be- 
tween the left and right sides) were: Hamiltonian 0.2%, momentum 3.3%. In the large-scale 
run these errors were smaller: Hamiltonian 0.05%, momentum 0.3%. We ran these runs 
also with a coarser grid. Going from a 32 3 grid to the 64 3 grid, the Hamiltonian error went 
down to about one-quarter, the momentum error to about one-half. 

VI. CONCLUSIONS 

We have created a working truly 3-dimensional numerical relativity code for studying 
inhomogeneous inflation. The code has been tested against perturbation theory with 
good results. On Cray Y-MP computers, a practical grid size to run the code with is 
around 48 3 to 64 3 . This is enough to handle an interesting amount of structure with 
reasonable accuracy (but just barely). Present memory limits would allow larger grids, 
but the increased computation time leads to much worse throughput. 

A variant of this code is being used[13] at present for interacting black hole studies, 
with initial data computed as described by Cook et al, [14] and with asymptotically flat 
(rather than periodic or reflective) boundary conditions. Sizes at least up to 128 3 are being 
implemented on a CM-5 computer[15]). We are working also to move the cosmology code 
to this machine where it will be run with inflaton and with other sources, to describe a 
number of 3-dimensional physical situations. 

One of the main motivations for the inflationary scenario was to explain the large- 
scale homogeneity and isotropy of the universe without requiring these properties as initial 
conditions[16]. An apparent shortcoming in the main body of work on inflation is then, 
that it has nevertheless been in the context of homogeneous space or small perturbations 
around it. We have now presented numerical simulations of inhomogeneous inflation, 
where the inhomogeneity has been truly 3-dimensional and nonperturbative. These runs 
had sufficient inflation to solve the standard cosmological conundra. In particular, they 
contain regions where evidence of the initial inhomogeneity would not be locally observable. 
We have thus demonstrated the viability of inflation with inhomogeneous initial conditions. 
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Figure captions 

FIG. 1. The inflaton field as a function of time in homogeneous chaotic inflation. 
The lines showing the code result and the exact result fall on top of each other and cannot 
be distinguished by eye. Inflation ends near t = 300tpi. 

FIG. 2. Same as Fig. 1, but for the scale factor a. Inflation ends after a linear 
expansion of 10 , or a volume expansion of 10 102 . 

FIG. 3. The amplitudes of the three perpendicular plane-wave perturbations 5(f>i, 
5(p2, and 8(ps in the inflaton field. The different line styles show code results with different 
grid resolutions, as well as the perturbation theory result, (a) Early part of the run, 
perturbations inside the horizon: The perturbations with the longer wavelengths have 
longer oscillation periods, but all are damped by the expansion at the same rate. In this 
case the code results are rather accurate, and the lines fall onto each other, (b) Late part 
of the run, perturbations outside the horizon: The vertical scale has been expanded by a 
factor of 25. The inset displays the last decade (from t = 10 3 tpi to t = 10 4 tpi) with a 
further expanded scale to show the post-inflation oscillations. 

FIG. 4. Same as Fig. 3, but for the perturbations Ai, A 2 , and A 3 in the metric, (a) 
Early part of the run. (b) The entire run. The perturbations exit the horizon between 
t = 2£pi and t = 3tp\ and the oscillatory behaviour changes to monotonic growth, which 
ceases as inflation ends. 

FIG. 5. Same as Fig. 3, but for the perturbations Bi, B 2 , and B3 in the metric. We 
only show the early part of the run, because in the late part B stays almost constant in 
time. 

FIG. 6. Regions of small and large values of <j> in the initial data. We express the 
contour levels as a percentage between the minimum (0%) and the maximum (100%) value 
in the slice, (a) The regions of small (p are enclosed by 20%, 40% and 50% contours. In 
the small-scale run these correspond to = 4.973mpi, 4> = 4.986mpi, and 4> = 4.993mpi. 
In the large-scale run the contours have <j> = 4.12mpi, <j> = 4.55mpi, and <j> = 4.77mpi. 
The minimum value, near the back corner enclosed by all three contours, was 4.959mpi 
in the small-scale run, and 3.69mpi in the large-scale run. (b) The regions of large <p 
are shown enclosed by 75% and 90% contours (5.010mpi and 5.020mpi for the small-scale 
run, 5.31mpi and 5.63mpi for the large-scale run). The maximum value was 5.026mpi in 
the small-scale run, and 5.85mpj in the large-scale run. In some of the following figures, 
quantities are plotted along a reference line. This line runs parallel to the y-axis, which in 
this figure goes to the right from the origin. The reference line goes through the initial <p 
minimum near the back corner and through the region with large <p near the right end of 
the y-axis. 

FIG. 7. Figures 7-14 show results from the small-scale run. This figure shows the 
inflaton field along the reference line at various time slices. These correspond to t/tpi = 5 
(initial time), 5.6, 6.4, 7.6, 9.2, 11.3, 14.3, 18.5, 24.1, 42.6, and 77.7. The Hubble time 
was initially lOtpj. At the last slice shown, where inflation has taken hold, the potential 
contributing 99.6% of total energy, the approximate Hubble time is 28tpi. 

FIG. 8. The inflaton field at t = 9.2tpi. We show regions of small <fi enclosed by 15% 
and 30% contours. 

FIG. 9. The local volume expansion Vh at t = 9.2£pi. (a) Re gions of small y/h are 
enclosed by 30%, 70%, and 83% contours, (b) Regions of large \fh are enclosed by 94% 
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and 98% contours. 

FIG. 10. Same as Fig. 7 but for \fh. 

FIG. 11. Trace of the extrinsic curvature at t = 9.2tpi. (a) Small values of \K\ are 
enclosed by 50%, 80%, and 88% contours, (b) Large values of \K\ are enclosed by the 98% 
contour. 

FIG. 12. Same as Fig. 7 but for K. 

FIG. 13. Scalar curvature at t = 9.2tpi. Regions of large R are enclosed by 25% and 
75% contours. 

FIG. 14. Same as Fig. 7 but for R. 

FIG. 15. Figures 15-20 show results from the large-scale run. This figure shows the 
inflaton field at cycle 512, or t = 166tpi, after about 6 Hubble times. Regions of small </> 
are enclosed by 15%, 34%, and 44% contours. 

FIG. 16. The inflaton field along the reference line at every 512th cycle. The initial 
slice is at the top, and the last slices at the bottom, as 4> moves down with time. The 
slices shown correspond to t/t ¥ \ = 5, 166, 482, 598, 852, 1140, 1471, 1861, 2337, 2947, 
3797, 5201, 9624, 6.02 x 10 4 , and 7.10 x 10 5 . Defining the end of inflation as the time when 
V((p) falls below 50% of the total energy, it happens at t = 1.2 x 10 4 tpp The run contains 
4 periods of post-inflation oscillations. The last two slices shown (close to each other) are 
from this era. 

FIG. 17. The local volume expansion \fh. (a) At cycle 64 or t = 14.7tpi. Regions of 
large \fh are enclosed by 60% and 80% contours, (b) At cycle 512 or t = 166tpp Regions 
of large \/h are enclosed by 5% and 50% contours. 

FIG. 18. Same as Fig. 16 but for \Jh. \fh moves up with time. 

FIG. 19. Same as Fig. 16 but for K. K moves up with time. The initial slice, where 
K is homogeneous, falls below the area plotted. 

FIG. 20. Scalar curvature at t = 14.7tpp (a) Regions of large positive R enclosed by 
15%, 30%, and 50% contours (R = 0.15 x 10- 3 m P1 ,0.55 x 10- 3 m P1 , 1.10 x 10- 3 m P1 ). (b) 
Regions of negative R enclosed by the 3.7% contour (R = —0.16 x 10 -3 ). 
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